home *** CD-ROM | disk | FTP | other *** search
/ Aminet 22 / Aminet 22 (1997)(GTI - Schatztruhe)[!][Dec 1997].iso / Aminet / mus / misc / MusicIn.lha / musicin / FastFT_float_asm.a < prev    next >
Text File  |  1997-09-23  |  15KB  |  509 lines
  1. ;--------------------------------------------------------------------------------------------
  2. ;--------        Fast Fourier Transformation for real-valued inputs           ---------------
  3. ;--           ©1997 Henryk "Buggs" Richter, tfa652@cks1.rz.uni-rostock.de                  --
  4. ;--                                                                                        --
  5. ;--        able to transform 1 real valued array of input data at highspeed                --
  6. ;--                                                                                        --
  7. ;--         Inputs & Outputs are single precision floating point variables                 --
  8. ;--                                                                                        --
  9. ;-- Requirements: 68020 + FPU                                                              --
  10. ;--                                                                                        --
  11. ;-- date: 17.09.97                                                                         --
  12. ;--------------------------------------------------------------------------------------------
  13. ;
  14.  
  15. Gamma        EQU        9        ;N = 2^Gamma (max input points this routine has got data arrays for)
  16.  
  17.                 section    blah,code
  18.  
  19.                   XDEF     @ASM_FastFT_main_loop
  20.                   XDEF     _ASM_FastFT_main_loop   ; for SAS-C 6.5
  21.                   XDEF     @ASM_FastFT_energy_phi
  22.                   XDEF     _ASM_FastFT_energy_phi  ; for SAS-C 6.5
  23.  
  24. ;-------------------------------------------------------------------------------------------
  25. ;          create Power spectrum + phase angle (based on Stéphane TAVENARD`s routine)
  26. ;
  27. ;Inputs: a0: float *x_real
  28. ;        a1: float *x_imag
  29. ;        a2: float *energy
  30. ;        a3: float *phi
  31. ;        d0: n
  32. ;
  33. _ASM_FastFT_energy_phi
  34. @ASM_FastFT_energy_phi
  35.                   movem.l  a2-a6,-(sp)
  36.                   
  37.                   lsr.w    #1,d0            ; process only n/2 + 1 data items (=HBLKSIZE)
  38. FFT_energy_phi_0
  39.                   fmove.s  (a0)+,fp0
  40.                   fmove.s  (a1)+,fp1
  41.                   fmove.x  fp0,fp4
  42.                   fmove.x  fp1,fp5
  43.                   fmul.x   fp4,fp4           ; (#2)
  44.                   fmul.x   fp5,fp5           ; (#2)
  45.                   fadd.x   fp4,fp5           ; fp5 = en = real^2 + imag^2
  46.                   fmove.s  fp5,(a2)+         ; *energy++ = en
  47.  
  48.                   ftst.x   fp5
  49.                   fbne.w   FFT_energy_phi_1  ; en <> 0
  50.                   fmove.s  #0,fp1            ; phi = 0
  51.                   bra.s    FFT_energy_phi_5
  52. FFT_energy_phi_1  ftst.x   fp0
  53.                   fbeq.w   FFT_energy_phi_3  ; real == 0
  54. FFT_energy_phi_2  fdiv.x   fp0,fp1           ; fp1 = imag / real (#2)
  55.                   fatan.x  fp1,fp1           ; fp1 = atan( imag / real )
  56.                   bra.s    FFT_energy_phi_5
  57. FFT_energy_phi_3
  58.                   fmove.s  #1.570796327,fp2  ; ¶/2 (90° in radiant)
  59.                   ftst.x   fp1
  60.                   fbgt.w   FFT_energy_phi_90
  61.                   fneg.x   fp2               ; -¶/2
  62. FFT_energy_phi_90                  
  63.                   fmove.x  fp2,fp1 
  64. FFT_energy_phi_5  fmove.s  fp1,(a3)+         ; *phi++ = phi
  65.                   dbra     d0,FFT_energy_phi_0
  66.                   
  67.                   movem.l  (sp)+,a2-a6
  68.                   rts
  69.  
  70. _ASM_FastFT_main_loop
  71. @ASM_FastFT_main_loop
  72.  
  73. ;-------------------------------------------------------------------------------------------
  74. ;Inputs: d0: power
  75. ;        a0: float *x_real    ;x_real[k]   with k=0,1,...,(1<<power)-1
  76. ;        a1: float *x_imag    ;x_imag[k]=0
  77. ;       (a1 as Input pointer irrelevant since only real-valued samples need to be transformed)
  78. ;        a2: float *SinCosTab
  79. ;        a3: float *SinCosTab2
  80. ;        (a5: WORD  *BitreverseTab)       ;at the moment this table kept local
  81. ;
  82. ;Outputs: x_real[k], x_imag[k] with k=0,1,...,(1<<power/2)-1
  83. ;         is the fourier transformed array[1<<power] only representing the
  84. ;         positive frequency range (the negative part is mirrored anyway
  85. ;         and wouldn`t contain extra information == redundant)
  86. ;
  87.         movem.l    d0-d7/a0-a6,-(sp)
  88.         subq.w    #1,d0            ;this FFT needs only half of
  89.                         ;the operations due to
  90.                         ;real valued input samples
  91.                         ;
  92.  
  93.         lea    BitreverseTab,a5    ;Bit reversing (Butterfly)
  94.  
  95.         cmp.w    #Gamma,d0        ;more data than we`ve got data arrays for ?
  96.         bhi    FFT_Fail
  97.         cmp.l    Last_Gamma,d0        ;need to create tables ?
  98.         beq    SkipTables
  99.         move.l    d0,Last_Gamma
  100.  
  101.         bsr    MakeBitreversetabs
  102. SkipTables
  103.         movem.l    d0/a0/a1/a3,-(sp)    ;need to be restored upon writing
  104.                         ;the N*2 transformed output out
  105.                         ;of the N points transformation
  106.  
  107.         bsr    MakeIntab        ;reorganize input tables for FFT:
  108.                         ;assumed, you wanna transform 1024
  109.                         ;successive points of the form:
  110.                         ;
  111.                         ;x(k) ; k=0,1,...,2N-1
  112.                         ;
  113.                         ;then you have to split the data
  114.                         ;in the following way:
  115.                         ;
  116.                         ;Realtab(n)=x(2n)   \ n=0,1,...,N-1
  117.                         ;Imgtab(n) =x(2n+1) /
  118.                         ;              
  119.         lea    Realtab,a0
  120.         lea    Imgtab,a1
  121. ;------------------------------- InitWerte ------------------------------------------------------
  122.         moveq    #0,d1
  123.         bset    d0,d1        ;1<<Gamma = N
  124.  
  125.         move.l    d1,d5
  126.         lsl.w    #2,d5
  127.         subq.l    #1,d5        ;N*4-1
  128.  
  129.         lsl    #1,d1        ;N2=N/2 (*4 because of longword-wide mem access)
  130.  
  131.         move.w    d0,d2        ;Gamma
  132.         subq.w    #1,d2        ;NU1=Gamma-1
  133.  
  134.         moveq    #0,d3        ;K=0
  135.  
  136.         move.l    a2,a3        ;Sine- and Cosine table for N points, bitreversed
  137.  
  138.         cmp        #2,d2        ;< 8 points to transform ?
  139.         blt        FFT_Standard
  140. ;-----------------------------------------------------------------------------------------------
  141. ;----------------------- 1st FFT run, optimized for sin=0 & cos=1 ------------------------------
  142. ;-----------------------------------------------------------------------------------------------
  143.         move.w    d1,a2        ;save d1, make DBF-Loop possible
  144.         lsr.w    #2,d1        ;d1 is originally 4 * N/2
  145.         subq.w    #1,d1        ;Loop from I to N2
  146.  
  147.         move.w    a2,d7        ;N2
  148.         add.w    d3,d7        ;K+N2
  149. FFTopt_Innerloop
  150.  
  151.         ;----- real part ---------
  152.  
  153.         fmove.s    (a0,d7.w),fp0    ;XReal(k+N2)
  154.  
  155.         fmove.s    (a1,d7.w),fp4    ;XImag(k+N2)
  156.         fmove.x    fp4,fp6
  157.  
  158.         fmove.s    (a0,d3.w),fp6    ;XReal(K)
  159.         fsub.x    fp0,fp6        ;minus Treal
  160.         fmove.s    fp6,(a0,d7.w)    ;XReal(K+N2)=XReal(k)-TReal
  161.  
  162.         fadd.x    fp0,fp6        ;-> =XReal(k)
  163.         fadd.x    fp0,fp6        ;Xreal(k)+TReal
  164.         fmove.s    fp6,(a0,d3.w)    ;write conjugated complex mirror result
  165.  
  166.         ;----- imaginary part ------
  167.  
  168.         fmove.s    (a1,d3.w),fp0    ;XImag(K)
  169.         fsub.x    fp4,fp0        ;minus TImag
  170.         fmove.s    fp0,(a1,d7.w)    ;XImag(K+N2)=XImag(k)-TImag
  171.  
  172.         fadd.x    fp4,fp0        ;-> =XImag(k)
  173.         fadd.x    fp4,fp0        ;XImag(k)+TImag
  174.         fmove.s    fp0,(a1,d3.w)    ;write conjugated complex mirror result
  175.  
  176.         addq.w    #4,d3        ;k=k+1
  177.         addq.w    #4,d7        ;k+N2=k+N2+1
  178.         dbf    d1,FFTopt_Innerloop
  179.  
  180.         move.w    a2,d1
  181.  
  182. ;------------------ one complete iteration done, next one -----------------------------------
  183.         moveq    #0,d3        ;K=0
  184.         lsr.w    #1,d1        ;N2=N2/2
  185.         subq    #1,d2
  186.  
  187.  
  188.  
  189. ;-----------------------------------------------------------------------------------------------
  190. ;---------- 2nd FFT run, optimized for sin=0 & cos=1 (and sin=1 & cos=0) --------------------
  191. ;-----------------------------------------------------------------------------------------------
  192.  
  193.         move.w    d1,a2        ;save d1, make DBF-Loop possible
  194.         lsr.w    #2,d1        ;d1 is originally 4 * N/2
  195.         subq.w    #1,d1        ;Loop from I to N2
  196.  
  197.         move.w    a2,d7        ;N2
  198.         add.w    d3,d7        ;K+N2
  199. FFTopt2_Innerloop
  200.  
  201.         ;----- real part ---------
  202.  
  203.         fmove.s    (a0,d7.w),fp0    ;XReal(k+N2)
  204.  
  205.         fmove.s    (a1,d7.w),fp4    ;XImag(k+N2)
  206.         fmove.x    fp4,fp6
  207.  
  208.         fmove.s    (a0,d3.w),fp6    ;XReal(K)
  209.         fsub.x    fp0,fp6        ;minus Treal
  210.         fmove.s    fp6,(a0,d7.w)    ;XReal(K+N2)=XReal(k)-TReal
  211.  
  212.         fadd.x    fp0,fp6        ;-> =XReal(k)
  213.         fadd.x    fp0,fp6        ;Xreal(k)+TReal
  214.         fmove.s    fp6,(a0,d3.w)    ;write conjugated complex mirror result
  215.  
  216.         ;----- imaginary part ------
  217.  
  218.         fmove.s    (a1,d3.w),fp0    ;XImag(K)
  219.         fsub.x    fp4,fp0        ;minus TImag
  220.         fmove.s    fp0,(a1,d7.w)    ;XImag(K+N2)=XImag(k)-TImag
  221.  
  222.         fadd.x    fp4,fp0        ;-> =XImag(k)
  223.         fadd.x    fp4,fp0        ;XImag(k)+TImag
  224.         fmove.s    fp0,(a1,d3.w)    ;write conjugated complex mirror result
  225.  
  226.         addq.w    #4,d3        ;k=k+1
  227.         addq.w    #4,d7        ;k+N2=k+N2+1
  228.         dbf    d1,FFTopt2_Innerloop
  229.  
  230.         move.w    a2,d1
  231.         add.w    d1,d3
  232.  
  233. ;------------------ 1st half of this iteration done ----------------------------------------
  234.  
  235.  
  236.     ;optimized for sin=1 , cos = 0
  237.         move.w    d1,a2        ;save d1, make DBF-Loop possible
  238.         lsr.w    #2,d1        ;d1 is originally 4 * N/2
  239.         subq.w    #1,d1        ;Loop from I to N2
  240.  
  241.         move.w    a2,d7        ;N2
  242.         add.w    d3,d7        ;K+N2
  243. FFTopt3_Innerloop
  244.         ;----- real part ---------
  245.  
  246.         fmove.s    (a0,d7.w),fp1    ;XReal(k+N2)
  247.  
  248.         fmove.s    (a1,d7.w),fp0    ;XImag(k+N2)
  249.  
  250.         fmove.s    (a0,d3.w),fp6    ;XReal(K)
  251.         fsub.x    fp0,fp6        ;minus Treal
  252.         fmove.s    fp6,(a0,d7.w)    ;XReal(K+N2)=XReal(k)-TReal
  253.  
  254.         fadd.x    fp0,fp6        ;-> =XReal(k)
  255.         fadd.x    fp0,fp6        ;Xreal(k)+TReal
  256.         fmove.s    fp6,(a0,d3.w)    ;write conjugated complex mirror result
  257.  
  258.         ;----- imaginary part ------
  259.  
  260.         fmove.s    (a1,d3.w),fp0    ;XImag(K)
  261.         fadd.x    fp1,fp0        ;minus TImag
  262.         fmove.s    fp0,(a1,d7.w)    ;XImag(K+N2)=XImag(k)-TImag
  263.  
  264.         fsub.x    fp1,fp0        ;-> =XImag(k)
  265.         fsub.x    fp1,fp0        ;XImag(k)+TImag
  266.         fmove.s    fp0,(a1,d3.w)    ;write conjugated complex mirror result
  267.  
  268.         addq.w    #4,d3        ;k=k+1
  269.         addq.w    #4,d7        ;k+N2=k+N2+1
  270.         dbf    d1,FFTopt3_Innerloop
  271.  
  272.         move.w    a2,d1
  273.  
  274. ;------------------ one complete iteration done, next one -----------------------------------
  275.  
  276.         moveq    #0,d3        ;K=0
  277.         lsr.w    #1,d1        ;N2=N2/2
  278.         subq    #1,d2
  279.  
  280.  
  281. FFT_Standard
  282. ;-----------------------------------------------------------------------------------------------
  283. ;--------------------------------- FFT routine main loop ---------------------------------------
  284. ;-----------------------------------------------------------------------------------------------
  285. FFT_Outerloop
  286. FFT_Innerloop2
  287.         move.w    d1,a2        ;save d1, make DBF-Loop possible
  288.         lsr.w    #2,d1        ;d1 is originally 4 * N/2
  289.         subq.w    #1,d1        ;Loop from I to N2
  290.  
  291.         move.w    a2,d7        ;N2
  292.         add.w    d3,d7        ;K+N2
  293.  
  294.         move.w    d3,d4        ;M = K % 2^NU1
  295.         lsr.w    d2,d4        ;
  296.         and.w    #~3,d4        ;
  297.  
  298.         ;---------------- Bit reversing P=IBR(M) already done --------------------
  299.         fmove.s    4(a3,d4.w*2),fp5    ;cos(M)
  300.         fmove.s    (a3,d4.w*2),fp7        ;sin(M)
  301.                     ;W^P = e^(2*PI*P/N)
  302.                     ;or   real=cos(2*PI*P/N)
  303.                     ;     img =sin(2*PI*P/N)
  304. FFT_Innerloop
  305.         ;----- real part ---------
  306.  
  307.         fmove.s    (a0,d7.w),fp0    ;XReal(k+N2)
  308.         fmove.x    fp0,fp1
  309.  
  310.         fmove.s    (a1,d7.w),fp4    ;XImag(k+N2)
  311.         fmove.x    fp4,fp6
  312.  
  313.         fmul.x    fp5,fp0        ;Xreal*cos
  314.         fmul.x    fp5,fp4        ;XImag*cos
  315.  
  316.         fmul.x    fp7,fp6        ;XImag*sin
  317.  
  318.         fadd.x    fp6,fp0        ;Treal= XReal*cos+Ximag*sin
  319.  
  320.         fmove.s    (a0,d3.w),fp6    ;XReal(K)
  321.         fsub.x    fp0,fp6        ;minus Treal
  322.         fmove.s    fp6,(a0,d7.w)    ;XReal(K+N2)=XReal(k)-TReal
  323.  
  324.         fadd.x    fp0,fp6        ;-> =XReal(k)
  325.         fadd.x    fp0,fp6        ;Xreal(k)+TReal
  326.         fmove.s    fp6,(a0,d3.w)    ;write conjugated complex mirror result
  327.  
  328.         ;----- imaginary part ------
  329.                     ;ximag*cos (see above) is in fp4
  330.         fmul.x    fp7,fp1        ;xreal(k+n2) * sin
  331.         fsub.x    fp1,fp4        ;timag = ximag*cos - xreal*sin
  332.  
  333.         fmove.s    (a1,d3.w),fp0    ;XImag(K)
  334.         fsub.x    fp4,fp0        ;minus TImag
  335.         fmove.s    fp0,(a1,d7.w)    ;XImag(K+N2)=XImag(k)-TImag
  336.  
  337.         fadd.x    fp4,fp0        ;-> =XImag(k)
  338.         fadd.x    fp4,fp0        ;XImag(k)+TImag
  339.         fmove.s    fp0,(a1,d3.w)    ;write conjugated complex mirror result
  340.  
  341.         addq.w    #4,d3        ;k=k+1
  342.         addq.w    #4,d7        ;k+N2=k+N2+1
  343.         dbf    d1,FFT_Innerloop
  344.  
  345.         move.w    a2,d1
  346.         add.w    d1,d3        ;k=k+N2
  347.         cmp.w    d5,d3        ;N*4-1 ?
  348.         blt.w    FFT_Innerloop2
  349.  
  350. ;------------------ one complete iteration done, next one -----------------------------------
  351.         moveq    #0,d3        ;K=0
  352.         lsr.w    #1,d1        ;N2=N2/2
  353.  
  354.         dbf    d2,FFT_Outerloop ;NU1=NU1-1 -> to NU1 = 0 proceeded, at <0 done
  355.  
  356. ;---------------------------- reorder the results (Butterfly) and ----------------------
  357. ;--- get the coefficients for the N*2 point FFT back from the coeffs of N point FFT ----
  358.  
  359.         movem.l    (sp),d0/a0/a1/a3    ;restore Data area for writing Output
  360.  
  361.         move.l    a0,a2
  362.         move.l    a1,a4
  363.  
  364.         lea    Realtab,a0
  365.         lea    Imgtab,a1
  366.  
  367.         moveq    #0,d1
  368.         bset    d0,d1
  369.         move.w    d1,d2            ;N
  370.         subq.w    #1,d1            ;N-1
  371.  
  372.         moveq    #1,d0
  373.         fmove.s    #0.5,fp6
  374.  
  375.         fmove.s    (a0),fp2        ;Xreal(0)
  376.         fmove.s    (a1),fp4        ;Ximag(0)
  377.         fadd.x    fp4,fp2            ;Xreal + Ximag
  378.         fmove.s    fp2,(a2)+        ;Gleichanteil
  379.         clr.l    (a4)+            ;Phasendrehung für Gleichspannung ist immer 0 -> Ximag = 0
  380. FFT_reorg
  381.         move.w    (a5,d0.w*2),d4        ;i=IBR(k)
  382.         move.w    (a5,d1.w*2),d5        ;i=IBR(k)
  383.  
  384.         fmove.s    (a0,d4.w*4),fp2        ;Xreal(n)
  385.         fmove.s    (a0,d5.w*4),fp3        ;Xreal(N-n)
  386.  
  387.         fmove.s    (a1,d4.w*4),fp4        ;Ximag(n)
  388.         fmove.s    (a1,d5.w*4),fp5        ;Ximag(N-n)
  389.  
  390.         fmul.x    fp6,fp2        ;divide all values by 2
  391.         fmul.x    fp6,fp3
  392.         fmul.x    fp6,fp4
  393.         fmul.x    fp6,fp5
  394.  
  395.         fadd.x    fp3,fp2            ;R(n)+R(N-n)
  396.         fmove.x    fp2,fp0
  397.         fsub.x    fp3,fp2
  398.         fsub.x    fp3,fp2            ;R(n)-R(N-n)
  399.         fmove.x    fp2,fp3            ;R(n)-R(N-n)
  400.  
  401.         fmul.s    (a3,d0.w*8),fp2        ;sin ¶n/N * ( R(n)-R(N-n) )
  402.         fmul.s    4(a3,d0.w*8),fp3    ;cos ¶n/N * ( R(n)-R(N-n) )
  403.         fsub.x    fp2,fp0
  404.  
  405.         fsub.x    fp1,fp1
  406.         fsub.x    fp3,fp1            ;-[ cos ¶n/N * ( R(n)-R(N-n) ) ]
  407.  
  408.         fsub.x    fp5,fp4            ;I(n)-I(N-n)
  409.         fadd.x    fp4,fp1            ;+I(n)-I(N-n)
  410.         fadd.x    fp5,fp4
  411.         fadd.x    fp5,fp4            ;I(n)+I(N-n)
  412.         fmove.x    fp4,fp5
  413.         
  414.         fmul.s    (a3,d0.w*8),fp4        ;sin ¶n/N * ( I(n)+I(N-n) )
  415.         fmul.s    4(a3,d0.w*8),fp5    ;cos ¶n/N * ( I(n)+I(N-n) )
  416.  
  417.         fsub.x    fp4,fp1
  418.         fmove.s    fp1,(a4)+        ;-[ sin ¶n/N * ( I(n)+I(N-n) ) ]
  419.  
  420.         fadd.x    fp5,fp0
  421.         fmove.s    fp0,(a2)+        ;+[ cos ¶n/N * ( I(n)+I(N-n) ) ]
  422.  
  423.         subq.w    #1,d1
  424.         addq.w    #1,d0
  425.         cmp.w    d2,d0
  426.         blt.w    FFT_reorg
  427.  
  428.         fmove.s    (a0),fp2        ;Xreal(0)
  429.         fmove.s    (a1),fp4        ;Ximag(0)
  430.  
  431.         movem.l    (sp)+,d0/a0/a1/a3    ;restore Data area for writing Output
  432.  
  433.         fsub.x    fp4,fp2            ;Xreal - Ximag
  434.         fsub.s    (a0),fp2
  435.  
  436.         fmove.s    fp2,(a2)+        ;Frequency N (Samplingrate / 2)
  437.         clr.l    (a4)+            ;Phase for fmax = 0 or 180° -> Imag always 0
  438.  
  439. FFT_Fail
  440.         movem.l    (sp)+,d0-d7/a0-a6
  441.         rts
  442. ;----------------------- Create the local Bitreverse table -----------------------------
  443. ;Input: d0: power (used by this FFT, input Power / 2)
  444. ;       a5: WORD *BitreverseTab
  445. ;
  446. MakeBitreversetabs:
  447.         movem.l    d0-d7/a0-a6,-(sp)
  448.         
  449.         moveq    #0,d5
  450.         bset    d0,d5
  451.         subq.w    #1,d5            ;N-1
  452.  
  453.         move    d0,d6            ;Anzahl der Bits
  454.         moveq    #0,d1
  455. BRT_Innerloop1
  456.         moveq    #0,d3
  457.         move    d1,d2
  458.         move    d6,d7
  459.         subq    #1,d7
  460. BRT_Innerloop2
  461.         roxr    #1,d2
  462.         roxl    #1,d3
  463.         dbf    d7,BRT_Innerloop2
  464.         move    d3,(a5)+
  465.  
  466.         addq    #1,d1
  467.         dbra    d5,BRT_Innerloop1
  468.  
  469.         movem.l    (sp)+,d0-d7/a0-a6
  470.         rts
  471.  
  472. ;----------------------- Create input data array for FFT -------------------------------
  473. ;Input: d0: power (used by this FFT, input Power / 2)
  474. ;       a0: float *x_real
  475. ;
  476. ;used Globals:
  477. ;    *Realtab
  478. ;    *Imgtab
  479. MakeIntab
  480.         movem.l    d0-d7/a0-a6,-(sp)
  481.  
  482.         lea    Realtab,a2
  483.         lea    Imgtab,a3
  484.  
  485.         moveq    #0,d1
  486.         bset    d0,d1
  487.         subq.w    #1,d1
  488. IN_tab
  489.         move.l    (a0)+,(a2)+    ;g(k)=x(2k)   k=0,1,...,N*2-1
  490.         move.l    (a0)+,(a3)+    ;h(k)=x(2k+1)
  491.         dbra    d1,IN_tab
  492.  
  493.         movem.l    (sp)+,d0-d7/a0-a6
  494.         rts
  495.  
  496.  
  497.         section    __MERGED,DATA
  498.  
  499. Last_Gamma        dc.l    0    ;check for nesessity of creating the tables
  500.     
  501.         section    __MERGED,BSS    ;rename to __MERGED for SAS
  502.  
  503. Realtab            ds.l    1<<Gamma ;N values (float)
  504. Imgtab            ds.l    1<<Gamma ;N values (float)
  505.  
  506. BitreverseTab        ds.w    1<<Gamma ;N values (WORD)
  507.  
  508.         END
  509.